Dynamics of drag and force distributions for projectile impact in a granular medium 
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Our experiments and molecular dynamics simulations on a projectile penetrating a two- 
dimensional granular medium reveal that the mean deceleration of the projectile is constant and 
proportional to the impact velocity. Thus, the time taken for a projectile to decelerate to a stop is 
independent of its impact velocity. The simulations show that the probability distribution function 
of forces on grains is time-independent during a projectile's deceleration in the medium. At all times 
the force distribution function decreases exponentially for large forces. 
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Craters on the earth and moon are similar to craters 
formed in laboratory experiments at much lower energies 
by using projectiles and explosives In labora- 

tory experiments at large impact energies, narrow jets 
have been observed to rise even higher than the initial 
height of the projectile 0,0- Recent experiments have 
determined how the shape, depth, and width of craters 
formed in granular media depend on the energy of the 
impact projectile [6j,|7j, but there is little known about 
the dynamics of a projectile during crater formation. 

We have studied the time evolution of projectile mo- 
tion. Our experiments and molecular dynamics simula- 
tions on a two-dimensional granular medium yield the 
time dependence of the drag force on projectiles. Sim- 
ulations for the same conditions are in accord with the 
experiment and also yield the time evolution of the forces 
on all of the particles; hence, we can study the time de- 
pendence of the force probability distribution function at 
different stages of the projectile motion. 

Our observations and simulations reveal three distinct 
regimes of the motion, as illustrated in Fig. ^ impact, 
where the projectile first hits the granular medium; pene- 
tration, where a transient crater forms and grains in front 
of the projectile are fluidized; collapse, where the pro- 
jectile has almost stopped and the deep transient crater 
collapses, forming a static crater that remains visible on 
the surface. 

Methods — In the experiment, a projectile of diame- 
ter D = 4.46 cm and mass 32.2 g was dropped into a 
bed of small particles (cylinders) contained between two 
glass plates with a separation 1.1 times the length of 
the cylinders. The initial projectile heights h (h < 80 
cm) correspond to impact velocities up to 400 cm/s. To 
reduce crystallization, two sizes of small particles were 
used: 12600 particles (84% of the total number) had di- 
ameter di = 0.456 cm (mass mi = 0.049 g) and 2400 
particles had diameter d 2 = 0.635 cm (mass m 2 = 0.097 
g) . To obtain a uniform granular bed with a reproducible 
area fraction before each drop of the projectile, the bed 
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FIG. 1: Snapshots of a projectile in the three distinct regimes 
of its motion in a bidisperse mixture of particles (cylinders). 
Experiment: The larger cylinders of the bidisperse mixture 
are colored black for visualization and are 40% larger in di- 
ameter than the grey cylinders. The dashed line shows the 
location of the surface after the collapse is complete. Simula- 
tion: the shading of each particle is proportional to the sum 
of the magnitudes of all the normal forces acting on that par- 
ticle; this renders visible the instantaneous force chains. The 
projectile is 9.8 times as large in diameter and 657 times as 
massive as the smallest particles. 



was fluidized with air flow that was slowly reduced to 
zero, yielding the same bed height (65c?i) and area frac- 
tion (81 ± 2%) for each projectile drop. The bed width 
was 225g?i. The position of the projectile, y(t), defined 
as the distance between the bottom of the projectile and 
the initial height of the bed, was determined with a high 
speed camera and a center of mass particle tracking al- 
gorithm ||. 

We modelled the system with a soft-core molecular dy- 
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namics (MD) simulation that used 15, 000 disks that had 
the same sizes and area fraction as the experiment. Any 
two disks (one of which can be the projectile) exert the 
following normal and tangential forces on one another: 
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F n = - [kS + m r ~f n \v n \9(v n )] n 
F[ s = min m r j s \v s \, (J.\F n \ s, 



(1) 
(2) 



where S is the length of overlap 0, Hoj , and v n and v s 
are the normal and tangential components of the surface 
velocity (n and s are unit vectors parallel to v n and v s ). 
The four parameters of the model were found empiri- 
cally for one impact velocity and the same parameters 
were used for all other simulations: k = 3.2 x 10 3 kg 
s -2 0, 0, 0, 0| is proportional to Young's modulus, 
7„ = 10 4 s _1 and j s — 8 x 10 3 s" 1 are viscoelastic con- 
stants, and = 0.28 is the static friction coefficient. rn r 
is the reduced mass (m,^ 1 — rri^ + nig 1 for two parti- 
cles A and B). The Heaviside function 9 in F n models 
an elastic-plastic interaction (e.g., see Fig. 8 of [l5j|h 
the use of the Heaviside function distinguishes our force 
model from previous soft-core MD simulations @, H3 |. 
Simulations with a more realistic form for F s yielded 
results not significantly different from our simple form, 
which is computationally more efficient. A comparison 
of the simulation output using time steps shorter than 1 
/is did not yield different results; a 1 /is time step was 
used in the results presented here. 

Results — The simulation results agree remarkably 
well with the laboratory observations, as Fig. [21 illus- 
trates. Both experiment and simulation reveal that the 
time taken for a projectile to slow to a stop in the granu- 
lar medium is independent of its velocity at impact. The 
large deceleration of the projectile at impact (see Fig. EI) 
is similar to that of a projectile incident on a liquid. How- 
ever, in contrast to the behavior of a projectile in a fluid 
17], in the granular medium there is a long penetration 
region in which the projectile's average acceleration is 
constant: y(t) is described by a parabola (Fig. EJa)), so 
v y (t) decreases linearly in time (Fig. [2J^b)) . Further, the 
acceleration is proportional to the impact velocity, as the 
inset in Fig.|2{b) illustrates: a y — avog, where the slope 
of the line gives a — 0.0064±0.0001 s/cm. Thus, the pro- 
jectile slows almost to a stop in a time t — 1/ag ~ 0.15 s, 
independent of vq. The projectile does not immediately 
come to a complete stop; rather it then moves very slowly 
downward over the next few seconds as the particles in 
the bed make small rearrangements in response to the 
collapse of the transient crater. 

The drag force on the projectile, while constant on the 
average, exhibits large fluctuations, which have a f~ 2 
spectrum (Fig. 0). 

The simulation determines all of the forces on each 
particle at every instance of time. Every force exerted by 
a particle on the projectile during a short portion of its 
travel is shown in Fig.0] At each point in the projectile's 




FIG. 2: (a) Position y(t) and (b) velocity v y (t) of the projec- 
tile as a function of time for different impact velocities, from 
both experiment (o) and simulation (solid lines). The two ver- 
tical dot-dashed lines give approximate boundaries between 
three regions: impact, where the projectile rapidly deceler- 
ates (cf. Fig. [3J; penetration, where the mean acceleration is 
constant, as illustrated by a dashed line fit in (a) of a parabola 
to the results from experiment and simulation for each vo; and 
collapse, where the projectile has almost stopped and the par- 
ticles above it are collapsing to fill the transient crater left by 
the penetration. The ordinate for (b) for each successive im- 
pact velocity vq < 363 cm/s is shifted by 30 cm/s for clarity. 
Inset: normalized acceleration of the projectile versus impact 
velocity from experiment (•) and simulation (<0). 



trajectory only a few particles exert a significant force on 
the projectile. Each peak in the magnitude of the force 
between an individual particle and the projectile in Fig.0] 
corresponds to a maximum force felt by the first particle 
in a force chain that extends downward. Each force 
chain consists of a string of particles in contact. The sum 
of the magnitudes of forces felt by each particle in this 
chain is much greater than the average for the particles in 
the bed, as can be seen in Fig. ^ (simulation), where dark 
chains of particles extend downward from the projectile 
into the particle bed. 

Results for the probability distribution P(F, t) of all 
normal forces between particles located in front of the 
projectile in a semicircular region of radius 1.5-D centered 
at the bottom-most point of the projectile are shown in 
Fig-El The distribution P{F, t) changes with time during 
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FIG. 3: The time series of the force on the projectile obtained 
from the simulation. The three regimes of motion are sepa- 
rated by dot-dashed lines. Inset: The power spectrum of the 
projectile acceleration during the penetration regime (0.02- 
0.15 s) for a projectile with initial velocity vo = 238 cm/s is 
described by f~ a with a = 2.1 ± 0.2. 
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FIG. 4: Vertical component of the force computed for every 
particle in contact with the projectile during part of the pen- 
etration regime (0.100 < t < 0.108 s in Figs. and 0. Each 
force grows, reaches a maximum (representing the inclusion 
of a particle in a particular force chain), and then decreases. 
Each type of line represents a particular particle; thus, the 
two arrows correspond to the same particle that appeared 
first at 12.344 cm and then reappeared at 12.501 cm. The 
projectile impact velocity was vq — 238 cm/s. The average of 
the total force on the projectile during this interval was 0.57 
N. 



impact but is time invariant during penetration: Fig. [S] 
shows the same distribution at times ti and £3, which are 
respectively early and late in the penetration regime. The 
presence of an inflection point F* in P(F, t) marks the 
beginning of exponential decay for large F. The cross- 
over to an exponential distribution at F* increases lin- 
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FIG. 5: Probability distribution of normal contact forces be- 
tween grains for a projectile with vo — 112 cm/s at the fol- 
lowing times: during impact (ti = 0.02s, •), early in the 
penetration regime (t<z = 0.05 s, ■), late in the penetration 
regime (tz = 0.12 s, +), and during collapse (ti = 0.20 s, A). 
The distribution decays exponentially for F > F* . The de- 
pendence of F* on the impact velocity is shown in the inset; 
the slope is 0.047 ±0.004 mN-s/cm. Each curve was obtained 
by averaging over 50 runs. 



early with v , as shown in the inset of Fig. EI After the 
projectile has almost stopped, the distribution is similar 
to that found in previous studies of equilibrium [19j and 
near equilibrium |20| force distributions . 

Discussion — Our experiments and simulations show 
that the mean drag force on a projectile dropped into 
a granular medium is constant during most of the pro- 
jectile's trajectory, and this drag force is proportional 
to the projectile's impact velocity. In our experiments 
inertia plays a major role. Interestingly, previous experi- 
ments with low constant velocities and negligible inertial 
effects also yielded a constant drag force in a granular 
medium [2l|. 

Since the deceleration of the projectile is proportional 
to the impact velocity (see inset Fig. |2Jb)), the projec- 
tile penetration depth is also proportional to the impact 
velocity. While our results are for a two-dimensional sys- 
tem, the linear dependence of the penetration depth on 
impact velocity has recently also been observed for pro- 
jectile impact in a three-dimensional granular medium 

The drag force on our projectile fluctuates strongly, as 
found also for cylinders dragged at small constant veloc- 
ities in experiments (v ~ 0.1 cm/s) |2l) and simulations 
(v ~ 2 cm/s) 0, 113, H3- The power spectrum of the force 
fluctuations has a /~ 2 dependence, as observed in mea- 
surements of fluctuations of the stress on a slowly sheared 
two-dimensional granular medium |2^j | and in measure- 
ments of the torque on a torsional pendulum in contact 
with a vibrofluidized granular bed [26|. The f~ 2 depen- 
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dence is explained by assuming random jumps in the drag 
force [25f . In our experiment these jumps originate from 
the variation of the forces exerted by the grains in contact 
with the projectile (Fig.^J. 

Finally, our simulations have yielded the normal con- 
tact forces for all particles in the bed. The distribution 
function for the forces on the particles in front of the 
projectile rapidly evolves immediately after the projec- 
tile makes contact with the bed, and then the distribution 
becomes stationary as the projectile penetrates the bed. 
This stationary distribution decays exponentially beyond 
an inflection point at F* whose value is linearly propor- 
tional to the impact velocity. This is the first determina- 
tion of the force distribution for a granular medium for 
an accelerating particle. During impact, our force dis- 
tribution is different from that measured for static beds 
|l9|. where the force distribution decayed exponentially 
at all times, as predicted by the q-model |2t| . 
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